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Abstract 

We present a SED model of dusty galaxies, in which the equation of radiative 
transfer is solved by assuming spherical symmetry. The temperature fluctuation of 
very small dust particles is calculated consistently with the radiative transfer. The 
adopted dust model consists of graphite and silicate grains and PAHs, whose relative 
fractions are determined for each MW, LMC and SMC type extinction curve. This 
model allows us to derive the intrinsic SEDs of stellar populations embedded in dusty 
ISM, which are very important indicators for the age of stellar populations. Therefore, 
the evolutionary phase of starburst galaxies which have frequently very dusty ISM 
can be investigated with this SED model. We show that the SEDs of Arp220 and M82 
can both be explained by the same single stellar population, despite the significant 
differences in the SEDs and the infrared luminosities. The apparent difference between 
their SEDs is mainly caused by the difference in the optical depth. In contrast, the 
SED of prototypical star-forming ERO, HR10, indicates that this galaxy is relatively 
old comparing to Arp220 and M82. It is found that, in the case of optically thin limit 
like elliptical galaxies, the optical depth cannot be inferred only from the SED, due 
to a degeneracy between the optical depth, galactic size, and the spatial distribution 
of dust; the latter two are important for estimating the average temperature of dust 
grains in elliptical galaxies. When the observed size of elliptical galaxies is adopted 
for the model geometry, SEDs can be used to constrain the spatial distribution of 
dust in elliptical galaxies. 

Key words: dust, extinction — galaxies: ISM — galaxies: starburst — galaxies: 
stellar content - radiative transfer 
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1. Introduction 



Age estimate of galaxies is one of the most important and essential tasks in modern 
astronomy, although difficult. Recent successive findings of high redshift galaxies allow age 
estimates of young galaxies to understand their formation processes. In this context, spec- 
trophotometry evolutionary models of galaxies have been widely applied to derive ages of the 
high redshift galaxies such as 53W002 (Windhorst et al. 1991), B2 0902+34 ( Lilly 1988; Bithell 
& Rees 1990; Windhorst et al. 1991; Eisenhardt & Dickinson 1992), 4C41.17 (Chambers & 
Chariot 1990; Chambers, Miley & van Breugel 1990; Windhorst et al. 1991; Graham et al. 
1994; Lacy & Rawlings 1996; Dey et al. 1997), 6C1232+39 (Eales et al. 1993), B2 0920+34 
(Chambers & Chariot 1990), and B2 1056+39 (Bithell & Rees 1990). However, the adopted 
photometric data for spectral energy distributions (SEDs) are mainly restricted to the UV- 
optical range in the rest frame, except for few galaxies. It turned out that in such cases the 
SED fitting is less helpful in order to determine uniquely the ages of these galaxies. This 
difficulty arises from the fact that the UV-optical SEDs of young stellar populations heavily 
obscured by dust are very similar to those of intrinsically old stellar populations; in other words, 
the UV-optical SEDs are degenerate in age and dust attenuation. In order to estimate the age 
of galaxies which contain a large amount of dust, it is necessary to analyse SEDs within a wide 
range of wavelengths including the dust emission (Takagi, Arimoto & Vansevicius 1999). 

The importance of detailed modelling of dust effects on SEDs of galaxies has been demon- 
strated recently (e.g. Mazzei, De Zotti & Xu 1994; Gordon, Calzetti & Witt 1997; Vansevicius, 
Arimoto & Kodaira 1997; Silva et al. 1998; Efstathiou, Rowan-Robinson & Siebenmorgen 
2000). However, the interpretation of the SEDs of dusty galaxies in the UV-NIR, which are 
the key to estimate the age of the stellar population (e.g. Takagi et al. 1999), is controversial. 
In the SED model by Silva et al. (1998), the SEDs of Arp220 and M82 in the UV-NIR are 
totally dominated by the relatively old stellar populations; i.e. unrelated to the very young 
stellar populations which are newly formed in recent starburst events. On the other hand, 
Efstathiou et al. (2000) show that the UV-NIR SED of M82 can be explained by the obscured 
young stellar populations. Although Silva et al. (1998) and Efstathiou et al. (2000) adopt 
the different photometric data with the different aperture size, this difference is not enough to 
explain the contradiction. Silva et al. (1998) divided the interstellar matter (ISM) of model 
galaxy into the star forming molecular clouds (MCs) and the diffuse medium. The formation 
process of stars is simulated, provided that stars are born only within the optically thick MCs 
and progressively escape from them. When the elapsing time from the onset of starburst is 
shorter than the time scale for stars to escape from MCs, all newly born stars are embedded in 
MCs and their light is completely hidden at the optical wavelengths, as they found for Arp220 
and M82. Therefore, in their model, the contribution from starburst stellar populations to 
UV-NIR SEDs is negligible for any aperture size of these galaxies. Note that Efstathiou et al. 
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(2000) derive the age of starburst by applying the SED model, while Silva et al. (1998) assume 
the age of starbursts to reproduce observed SEDs. 

The following observational results provide important clues to the origin of the UV- 
NIR light of starbursts: 1) Rakos, Maindl & Schombert (1996) showed that the UV-optical 
SEDs of starburst galaxies are dominated by strongly reddened, young stellar populations; 2) 
Schmitt et al. (1997) demonstrated that UV-NIR SEDs of starburst galaxies are different from 
those of normal galaxies; 3) Meurer et al. (1997) showed a tight correlation between the UV 
spectral slope and the ratio of FIR flux to UV flux; 4) Goldader et al. (1995) showed that the 
2 /im continua of the ten ultraluminous infrared galaxies are dominated by red supergiants in 
starburst stellar populations, by using the CO indices of these galaxies. All these observations 
indicate a significant contribution of the starburst population to the UV-NIR SEDs of starburst 
galaxies. Therefore, in this paper, the UV-NIR SEDs consistently coupled with FIR-submm 
SEDs are used to estimate the age of starburst galaxies. In this way, it should be possible to 
reduce the uncertainties of estimated age considerably. 

The SED of starburst galaxies has sometimes been modelled by using the empirical 
relation between IRAS colours and total IR luminosities (e.g. Devriendt, Guiderdoni & Sadat 
1999; Totani & Takeuchi 2002). We would like to point out a vital importance of solving the 
radiative transfer problem in SED modelling, since it is the only method to employ the relations 
between the SED features in the whole wavelength range and physical parameters of the galaxies 
and dust consistently. Therefore, some additional spectral information is very helpful in order 
to constrain physical parameter space allowed for fitting, e.g. the silicate absorption feature 
at 9.7 /im arising due to the self-absorption by dust, can be used to constrain the type of 
extinction curve in our model. 

Here, we present a SED model of dusty galaxies by solving the equation of radiative 
transfer, in which stars are diffusely distributed within a spherical region. This assumption for 
geometry is clearly different from those by Efstathiou et al. (2000) in which a starburst region is 
modelled as an ensemble of star-forming molecular clouds. The Hubble Space Telescope (HST) 
revealed that most of nearby starburst galaxies actually have clusters of newly formed massive 
stars (O'Connell, Gallagher & Hunter 1994; Meurer et al. 1995; Tacconi-Garman, Sternberg & 
Eckart 1996), so-called super star clusters (SSCs). However, Meurer et al. (1995) reported that 
the SSCs contribute only 20% of the UV luminosity on average. Thus, UV light from starburst 
region is dominated by the light from diffusely distributed stars as modelled here, rather than 
that from star clusters. Also, note that Shioya, Taniguchi & Trentham (2001) suggest that in 
Arp220 the contribution from SSCs to the bolometric luminosity is at most 20%. 

Even if the observed SED of a galaxy is reproduced by the SED models, it is important 
to verify the uniqueness of the solution. Generally, a sophisticated handling of the SED model 
of dusty galaxies requires a number of parameters. In such a case, all the parameters may not 
be determined from the SED fitting alone; thus, the solution may not be unique. Therefore, it 
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is important to investigate the dependence of SEDs on each parameter. Once these test results 
are obtained, it becomes clear which parameters can be determined from the SED fitting. 

In this paper, we first intend to study the basic properties of SEDs of dusty galaxies, 
which are helpful to have an insight on what physical parameters can be constrained with the 
specific properties of SEDs. Then, we apply our SED model to dusty galaxies, a representative of 
the most active starburst (Arp220), a nearby less luminous starburst (M82), a distant extremely 
red star- forming galaxy (HR10), and an elliptical galaxy (NGC2768) as a study for optically thin 
case. Specifically, we focus on the starburst age of dusty starbursts derived from UV-submm 
SEDs. 

This paper is organized as follows. The model description and the SED analysis of 
parameter sensitivity are given in §2 and §3, respectively. Then, we apply our model to the 
observed SEDs in §4. Discussion and conclusions are given in §5 and §6, respectively. The 
cosmological parameters adopted throughout this paper are H = 75 km s" 1 Mpc" 1 , = 1 and 
A = 0. 

2. Model 

We have developed a code which solves the equation of radiative transfer by assuming 
spherical symmetry for arbitrary radial distributions of stars and dust. Isotropic multiple 
scattering is assumed and the self-absorption of re-emitted energy from dust is fully taken 
into account. The adopted dust model takes into account graphite and silicate grains, and 
polycyclic aromatic hydrocarbons (PAHs). We successfully reproduce the typical extinction 
curves in Milky Way (MW), Large and Small Magellanic Cloud (LMC and SMC) as well as 
the spectra of cirrus emission in the MW, which are modelled by considering the temperature 
fluctuation of very small grains and PAHs. As an input to our modelling, the intrinsic stellar 
SEDs of galaxies are taken from Kodama & Arimoto (1997). We assume that the SEDs of 
galaxies are dominated by the starburst region or by the diffusely distributed stars and dust: 
the former case is for starburst galaxies and the latter one for elliptical galaxies. In both cases 
the adopted geometry is defined by the radial distribution of stars and dust with different scale 
lengths. 

2.1. Stellar Spectral Energy Distribution of Galaxies 

An evolutionary synthesis code developed by Kodama & Arimoto (1997), which gives 
consistent chemical and spectral evolution of a galaxy, is used for generating the intrinsic stellar 
SEDs. The basic structure of the code follows Arimoto & Yoshii (1986) prescription, i.e. the 
effects of stellar metallicity are explicitly taken into account in calculating galaxy spectra. The 
characteristics of this model are as follows. 

• New stellar evolutionary tracks are incorporated comprehensively (VandenBerg et al. 1983; 
Bressan, Chiosi, & Fagotto 1994; Fagotto et al. 1994a, 1994b) and late stellar evolutionary 
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stages (horizontal branch - Yi 1996; asymptotic giant branch [AGB] - Vassiliadis & Wood 
1993; post-AGB and white dwarf - Vassiliadis & Wood 1994) are fully included 

• For the stellar spectra, Kurucz's (1992) stellar flux library is used for the metallicity 
0.0001 < Z < 0.05, the effective temperature 4000 K < T cS < 50000 K, and the surface 
gravity 0.0 < logg < 5.0. The library covers the wavelength range from 9 nm to 3.2 /im 
with the spectral resolution AA = 1 nm - 2 nm for A < 1 /im and AA = 5 nm - 10 nm for 
A > 1 um. 

• The spectra of cool stars with T e g < 4000 K are taken from Pickles (1985) and Gunn & 
Stryker (1983). 

Throughout this paper, we adopt the model of simple stellar populations (SSPs) with 
solar metallicity as the intrinsic stellar SED model, otherwise noted explicitly. We adopt the 
Salpeter initial mass function (IMF) with stellar masses ranging from 0.1 M Q to 60 M . Some 
examples of SEDs of SSPs are shown in Figure 1. In general, the intrinsic SEDs of stellar 
populations are mainly determined by both the age and the star formation history. Therefore, 
in order to derive the definite age of stellar populations from the SED, it is necessary to 
determine the time scale of star formation, which can be expressed in terms of the dynamical 
time, sound-crossing time, and cooling time. In this study, we mainly focus on starburst 
galaxies, in which the variation of stellar populations is less significant, comparing to normal 
spiral galaxies. Even in such an extreme case, it is difficult to know the characteristic time scale 
of each starburst event. Therefore, we choose the simplest stellar population as an indicator 
of the age of starbursts, rather than try to determine precise ages of starbursts from intrinsic 
stellar SEDs. 

Gas emission is not taken into account in our model which is designed to calculate 
the continuum emission of galaxies. A modification of SEDs due to the gas emission is not 
significant, unless a galaxy is younger than 10 Myr for SSPs (Leitherer & Heckman 1995; Fioc 
& Rocca-Volmerange 1997). As we will show later, our model can reproduce the observed SEDs 
of dusty galaxies with the age of SSP older than 10 Myr. Therefore, we believe that the effect 
of gas emission does not affect our conclusions. 

2.2. Dust model 

2.2.1. Optical Properties 

The optical properties of graphite and silicate are taken from Laor & Draine (1993). We 
adopt a " smoothed astronomical silicate" which is obtained by removing the absorption feature 
at 1/A = 6.5 urn -1 from the imaginary part of the original dielectric function (Draine & Lee 
1984; Laor & Draine 1993; Weingartner & Draine 2001). 

Optical properties of PAH are adopted according to Leger, d'Hendecourt & Defourneau 
(1989), and Desert, Boulanger & Puget (1990). The PAH features are assumed to have the 
Lorentz profile, which successfully reproduced the spectra of the reflection nebula NGC7023 
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and the Ophiuchus molecular cloud observed by ISO (Boulanger et al. 1998). The Lorentz 
function is given by: 



where f(y) is the normalized Lorentz profile, and A is the integrated cross-section, and uq and 
T give the frequency and the width of the band, respectively. The adopted parameters for the 
PAH features at A = 3.3, 6.2, 7.7, 8.6, 11.3 and 12.7 /im are given in Table 1. The number of H 
atoms attached to a PAH molecule is given by Nh = £#(6 Ac) 1 / 2 , where xh (< 1 by definition) 
is the degree of dehydrogenation, and Nc is the number of carbon atoms in one PAH molecule 
(Omont 1986). Since the binding energy of C-H bond is weaker than that of C-C bond, the H 
atoms tend to be removed in the strong radiation field. This suggests that xh could vary with 
the strength of radiation field. However, recent observations showed that relative intensities of 
PAH features are independent of the strength of radiation field (Lemke et al. 1998; Chan et al. 
2001). Thus, it is too early to introduce the relation between xh and the strength of radiation 
field, instead we assume xh = 0.4 according to Desert et al. (1990). 

2.2.2. Size Distribution 

The spectra of diffuse high latitude clouds (so-called Galactic cirrus), obtained by the 
IRAS all-sky survey (Low et al. 1984), show excesses of emission at 12 and 25 /mi over those 
expected from dust grains radiating at the equilibrium temperature. Draine & Anderson (1985) 
suggested that the excessive emission at MIR could be radiated by grains small enough to 
undergo the temperature fluctuation in the local interstellar radiation field, due to their small 
heat capacity. Such very small grains (VSGs), typically smaller than 0.01 /mi, are transiently 
heated up to several hundreds of Kelvin by the interstellar radiation field and emit in the MIR. 
Hereafter, we define both graphite and silicate grains with the radii a < 0.01 jum and a > 0.01 /mi 
as VSGs and big grains (BGs), respectively. 

For BGs, we adopt a power law size distribution derived by Mathis, Rumpl & Nordsiek 
(1977, hereafter MRN). For VSGs, the same size distribution is adopted, but with the different 
index (Draine & Anderson 1985). Introducing the radius = 0.01 /im, at which a break 
between the two power laws takes place, the size distribution is given by: 



where dnk is the number density of grains of type k (k = g for graphite and k = s for silicate) 
in the interval [a,a + da], uh is the number density of H nuclei, and C& is the coefficient which 
determines the ratio of dust to hydrogen. For both graphite and silicate grains, we adopt 
(3 = —3.5, a min = 0.001 /mi, a b = 0.010 /im and a max = 0.25 /im. The indices of VSG size 
distribution 7& are determined by the fits to the MW extinction curve and the spectrum of 
Galactic cirrus. 




(1) 




nnCka^a/ab^da, 
nnCka^da, 



for a min < a < a b (VSGs) 
for a b <a< a max (BGs) 



(2) 
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For PAHs, we write the size distribution as a power law of Nc'- 
drip = n H C P N^ p dN c , for N C1 < N c < N C2 , 



(3) 



where Nci and Nc2 are the number of C atoms in the smallest and largest PAH, respectively. 
We use the subscription k = P for PAHs. Since the small PAHs with Nc < 20 will evaporate 
even in the weak radiation field like the local interstellar radiation field (Omont 1986), Nci = 20 
is adopted in this work. On the other hand, Nc2 is less constrained. Nc2 should be less than 
~ 1000 so that a single photon could heat up PAHs to the temperature at which the emission 
peak is located at MIR (Puget & Leger 1989). We assume that the smallest graphite grain and 
the largest PAH contain the same number of C atoms (see also Dwek et al. 1997). According 
to the relation between a and Nc'- 



we get Nc2 = 465 for a m i n = 0.001 /im. It is assumed that 7p is equal to the index of size 
distribution of graphite VSGs, which is written as a function of Nc- Using above relations 
between a and Nc, 7p is given by (j g — 2)/3. 

2.2.3. Temperature Fluctuation 

When the dust grains attain temperature equilibrium with a surrounding radiation field, 
J\, by radiating the same amount of energy as the absorbed one, the equilibrium temperature, 
T eq , is determined by: 



where cr^ k (a) is the absorption cross-section of dust grains of type k and size a, and B\ is the 
Planck function. 

When dust particles are small enough, grains cannot attain the temperature equilibrium 
with the interstellar radiation field. In this case, we have to calculate the temperature distribu- 
tion function dP/dT, which is the probability for finding a particle in the interval [T,T + dT], 
in order to predict the emission spectra. Methods to calculate dP/dT have been presented 
by several authors (Draine & Anderson 1985; Desert, Boulanger & Shore 1986; Dwek 1986; 
Guhathakurta & Draine 1989). We follow the method outlined by Guhathakurta & Draine 
(1989), in which the transition matrix is defined as the probability per unit time of a particle 
making a transition in the enthalpy space by the emission or absorption, and the solution is 
derived by inverting this transition matrix. The heat capacities of graphite and silicate are also 
taken from Guhathakurta & Draine (1989). For PAHs, the same heat capacity is used as that 
of graphite. Although the different expressions are given for its heat capacity (e.g. Leger et al. 
1989; Dwek et al. 1997), alternative choice has no significant influence on the emission spectra 
(Siebenmorgen & Kriigel 1992). By using dP/dT, the emissivity per dust particle with the size 
a and the constituent k is given by: 




1/2 



for graphite grains 
for PAHs 



(4) 




(5) 
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B x (T)-±(a,T)dT, 



(6) 



where the sublimation temperature T su ^k = 800 K is assumed for all the dust constituents 
according to Carico et al. (1988) who showed that the emission from hot-dust particles with 
the temperature of ~ 800 K contributes to the NIR luminosity of IRAS bright galaxies. This 
assumption is safe because, unless T su ^k is higher than 800 K and the starburst is younger than 
10 Myr, the resulting SEDs depend little on T subjk . If we adopt T subik =1500 K (c.f. Kobayashi 
et al. 1993), a little excess of emission (< 10%) can be seen in the NIR wavelengths only for 
the optically thick cases in which dust emission dominates the NIR emission, and therefore this 
effect is negligible to the results of SED fitting. In case of the temperature equilibrium, dP/ dT 
is given by the delta function 5{T — T eq ). 

2.2.4. Dust in the MW, LMC and SMC 

Figures 2 and 3 show the model fits to the extinction curve and the spectrum of cir- 
rus emission in the MW, respectively. The extinction curve is given by the cross-section per 
hydrogen atom: 



where c\ k {a) is the extinction-cross section of the type k and size a dust grains. The adopted 
interstellar radiation field is taken from Mathis, Mezger & Panagia (1983). The contribution 
of PAHs to the spectrum and the extinction curve is significant; they dominate the spectrum 
at A < 12 /im, and are responsible for the nonlinear extinction rise at FUV. The predicted 
extinction at FUV is a factor of two higher than the average extinction curve in the MW. The 
dust model presented by Dwek et al. (1997), in which the same optical properties are adopted, 
also overestimates the FUV extinction curve. Unfortunately, as far as the visible- UV absorption 
of PAHs in the gas phase is concerned, very few studies are available especially for the large 
PAHs (Verstraete & Leger 1992). We think, this discrepancy should be solved by using the 
more reliable optical properties of PAHs than those available at present. 

Since the cirrus spectra of the LMC and SMC are not available at present, only the 
extinction curves are used to derive relative grain fractions, provided that fraction ratios of 
carbonaceous particles (graphite BGs, VSGs and PAHs) are the same as those of the MW (see 
also Pei 1992). The model fits to the extinction curve in the LMC and SMC are also shown 
in Figure 2. For normalisation of the extinction curves, the adopted ratios of hydrogen column 
density to colour excess E(B - V) are 5.8 10 21 , 2.0 10 22 , and 4.6 1 22 for the MW, LMC, 
and SMC, respectively (Mathis 1990). In Table 2, we give the derived mass ratios of dust to 
hydrogen and the relative mass fractions of various constituents of dust. Note that the fraction 
of VSGs of type k, f k SG and that of BGs, f k G , are related by: 




(7) 



fZ SG 

fBG 
Jk 




(8) 



Ik + 4 ai + a i - 4 + 4 ' 
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It is difficult to determine the amount of silicate VSGs, because they have little influence on 
both of the extinction curve and the spectrum of cirrus emission. Therefore, we simply extend 
downward the lower end of MRN distribution from a& to a m i n for silicate, i.e. 7 S = /3. For 
graphite VSGs, we derive 7 9 = —3.75. 

2.3. Radiative Transfer 

2.3.1. Geometry 

The density distributions, p(r) , of stars and dust are assumed to be given by King's law: 
p(r) = ~ TT' 0) 

where po is the density at the centre of galaxy and r c is a core radius. We fix a geometry of 
stars and change the dust distribution by choosing a value of parameter p = t c d/t c ,Si where 
r Cj 5 and r c ,D are the core radius of stars and dust, respectively. 

The second parameter that defines the geometry is a tidal radius r t which determines 
Po with a given total mass. When the total mass of stars is fixed, the temperature of dust 
is controlled by the tidal radius of stars r^s-, since the energy density of the radiation field is 
approximately proportional to the stellar density. In order to keep a source function (see below) 
constant, we adopt a mass-radius relation: 

r -^ 10R (JkT lkpc1, (10) 

where M* is the total stellar mass and R is a scaling factor, and the normalisation is given 
for the scale of typical galaxy. Note that, for constant value of R, the relative shape of SED 
is conserved for different M* due to this relation, when the optical depth is kept constant. 
We assume that the tidal radius of dust r t D is the same as r^s, i.e. r t D =r t s (=r t hereafter), 
and there is no dust at r > r t . We adopt the concentration parameter c = log(r t /r C] s) = 2.2, 
which is the typical value for normal elliptical galaxies (Combes et al. 1995). This value is 
also adopted for starburst galaxies in which the real distribution of stars inside the dominant 
starburst regions is difficult to derive due to the effect of dust attenuation. However, the effects 
of c on the SED are marginal, and detectable only at MIR region, and therefore have negligible 
impact on the values of parameters derived in §4. 

The dust distributions pu(r) for p — 1, 10, 100, and 1000 are shown in Figure 4. Note 
that, for larger values of p, r t can be smaller than r C) o, due to the definition given above. In 
the case of p = 1000, dust distributes nearly uniformly and surrounds centrally concentrated 
stars, while in the case of p = 1, stars and dust are distributed identically. The former case 
is similar to the shell geometry suggested for starburst galaxies by Gordon, Calzetti & Witt 
(1997). However, note that the shell geometry gives SEDs different from our case if p = 1000 
and Ty 3> 1, since the shell geometry results in the spectral cut-off around NIR (see the results 
by Rowan- Robinson & Efstathiou 1993). 

9 



A mass of dust, Mp, can be derived as follows. The optical depth at V band is defined 

by: 

m 

TV = / n H (r)a v dr, (11) 
Jo 

where nn{r) is the number density of hydrogen and o~y is the extinction cross-section of dust 
per hydrogen at V band. A total number (or total mass) of hydrogen is derived by using the 
distribution function defined by 77, once ry is specified. The mass ratio of dust to hydrogen 
is determined by the dust model. Since the column density is proportional to r t ~ 2 for a given 
mass, Md is given by: 

M ^ M ^fe) 2 - < 12 > 

where both M and r are the normalisation factor. By using the dust model described in 
§2.2 for the MW, LMC and SMC, M is found to be 5.5 10 7 M o , 9.4 1O 7 M and 1.1 1O 8 M , 
respectively. The second normalisation factor Tq depends on the geometry and is equal to 1.0, 
1.78, 33.86, and 1753 for 77 = 1000, 100, 10, and 1, respectively. 

2.3.2. Method of Calculation 

Solving the equation of radiative transfer for an arbitrary radial distribution of stars and 
dust, we adopt the "method of rays" (Band & Grindlay 1985). The schematic picture of this 
method is shown in Figure 5. 

It is convenient to divide the specific intensity, I\, into three components, 

W<» + /?> + /?>, (13) 

where the superscripts (1), (2) and (3) represent direct light from stars, scattered light by 
dust grains, and thermal emission from dust grains, respectively (Rowan- Robinson 1980). In 
the case that the optical properties of dust grains are independent of their temperature, each 

component of intensity can be calculated independently in the order of the superscript. The 

(2) (3) 

details of iterative calculation of I x and I x are discussed in the following subsections. 
The differential equation for the specific intensity, I x , along a ray is given by: 

= -n B (r)v? [h(b,s) -S x (r)}. (14) 

The impact parameter to the centre, b, and path length along the ray, s, are related to the 
radial distance, r, according to an expression: 



r 2 = b 2 + ^r 2 t -V-sf. (15) 
The source function for I x is given by: 

S x l \r)= A g , (16) 

47m H (r)af' 

where e* x (r) is the emissivity of stellar component at wavelength A. 
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Spherical geometry of the model requires only radial grid to be specified. The grid is 
defined by: 

r ' = (£r)' r * (17) 

where N is the total number of radial shells and d is a parameter that depends on the dust 
distribution. The required N depends on the adopted ry. In all analysed cases, the grid with 
d = 3 and iV = 160 gave results with satisfying accuracy; i.e. the difference between the input 
and output energy is less than 5% (see §2.3.4). If the required accuracy is violated the finer 
grid should be redefined. 

For each bj (= r^), 1^' is calculated at the crossing points of the ray and shells' bound- 
aries. The grid of path length, Sj, is defined at these points. Then, the equation for j'-th ray is 
written as: 

7i 1) (6„ Sm ) =§ i + e- A ^(7l 1) (6 i , Si ) - (18) 

with the boundary condition 7^ (fey,si) = 0, where Arj = (s i+1 — s^n^er^, and the bar above 
variables indicates an average of i-th and (i + l)-th values. Once I\(bj,Si) is calculated for all j 
(1 < j < N) and % (1 < % < 2(N — j) + 1), one can derive the mean intensity in case of spherical 
geometry: 

4 1) W = ^£ 1 7i 1) (r, yU )d / i, (19) 

or 

J ( x\r j/ ) = l^(I^(b J ,s J+ ) + I?(b,,s^))A H , (20) 

Z 3=1 

where A/ij is the step of cosine grid, fij, which is defined by the angle between the vector 
pointing the direction of the j-th ray through the radial point ry (y in Figure 5) and the radius 
vector; sj + and Sj- are the path lengths of the ray defined by bj, from both sides of the outer 
edge to the shell ry, which means j+ = N + (j' — 1) — 2(j — 1) and j— = N — (j' — 1); the bar 
above 7a means average of the j-th and (j + l)-th values. 

The components, 1^ and are calculated in a similar manner by using the source 
function characteristic for each component. 

2. 3. 3. Scattered Light 

The source function of n-times scattered light is: 

S?l(r,n) = ^-J »(m,n)^ti(r, m)dO, (21) 

where 7{ 2 ^_ x is the intensity of (n — l)-times scattered light, g(m,n) is the angular phase 
function for coherent scattering from direction m to n, and u\ is a dust albedo (e.g. Mihalas 
1978). Isotropic scattering is assumed for simplicity, that means g(m,n) = 1. The total 
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intensity of scattered light, 1^ is given by J2 n I\n- Since 1$/ can be regarded as the 0-times 
scattered light, the source function of once scattered light in the case of spherical geometry is: 

Sg(r) = £/ Jg(r,n)<n = f £' J?>(r,^ (22) 

Using this definition of the source functions, the intensities of the higher terms of scattered 
light are calculated by using the method given in the previous section. When the number of 
scattering terms considered is not sufficient, this results in a loss of the output energy. We find 
that the losses are less than 0.3% for the MW extinction curve when up to six terms are taken 
into account. 

In Figures 6 and 7, the mean intensity ratios of n to n — 1 times scattered light as a 
function of r are shown for different t] and ry in the case of uj = 0.6. Generally, the mean 
intensity ratio approaches to a flat distribution for the higher order of scattered light. For large 
Ty, this mean intensity ratio is hard to become totally constant from r = to r t , due to the slow 
diffusion of photons. Figure 8 shows the fraction of the unabsorbed stellar light (direct light) 
and the emerging scattered light both normalized with respect to the total input luminosity 
for the uniform distribution of stars and dust (Table 3A in Witt, Thronson & Capuano 1992). 
In our calculation, both values of the direct and scattered light are systematically smaller in 
comparison with the results obtained by Witt et al. (1992) who solved the radiative transfer 
problem by using the Monte Carlo method. The analytic solution for the direct light in an 
optically thick case is given by Band & Grindlay (1985) and is plotted together. Our results 
are consistent with the analytic solution. In Figure 9, luminosity ratios of scattered to direct 
light are compared. The results of Witt et al. (1992) are plotted by circles, while ours are 
represented by solid lines. Each solid line corresponds to a particular number of scattering 
terms used for calculation of the scattered light. Over the whole range of Ty, the results of 
Witt et al. (1992) are consistent with ours when the first five scattering terms are taken into 
account. 

2.3.4. Self-absorption 

In order to take into account a self-absorption of dust emission, the radial distribution 
of dP/dT, is iteratively calculated for every constituent and every discrete grain size until 
the energy conservation is fulfilled. The first guess of dP/dT, which is denoted as dP/dTi, is 
done using the mean intensity calculated from only two components, /j; and 1^ . With the 
assumption that the scattering of dust emission is negligible, the source function for the n-th 

(3) 

guess of the thermal component of intensity, ij^, is given by: 

Sg(r) = £ / < k (a)B x (T)^(a,T,r)dT^da. (23) 

n H {r)af ^ J Jo dT n da 

Calculation of the 1^ is performed in the same way as described in §2.3.2. The sum of 
intensities, 1^ + 1^ + I\l is used to derive dP/dT n+ i. The mean temperature of dust rises step 
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by step during the iterative procedure. We continue the iteration until the energy conservation 
is fulfilled. Even in the worst cases (ry ~40), the grid can be adjusted to keep the "numerical" 
losses of the total energy below 5%. This numerical loss is less than 3% in the nominal case. 

The changes of the spatial distribution of equilibrium temperature of BGs, T eq (r), and 
the SED due to the self-absorption are shown in Figure 10, panels a and b, respectively. Heating 
dust by the thermal emission of dust itself, i.e. the self-absorption, has considerable effect on 
the MIR emission if ry is larger than ~ 10. The effect of self-absorption is most prominent in a 
case of the SMC extinction curve since it has the largest absorption cross section at MIR due 
to high fraction of silicate grains. 

Ivezic et al. (1997) have defined a set of benchmark problems and gave exact solutions 
suitable for verification of the radiative transfer codes dealing with dusty environments. The 
benchmark problems are solved with a central point source embedded in a spherically symmetric 
dust envelope with an inner cavity being free of dust. The density variation with r is assumed 
to be a power law, p(r) oc r~ p . We show the results of comparison with p = in Figure 11. 
The consistency is very good although, in the extreme case of T\^ m = 1000, the temperature 
of dust is slightly higher than the benchmark solution. By the definition of the benchmark 
problem, the radius of inner cavity, r 1; has to be adjusted during the iteration to keep the dust 
temperature at the inner edge of dust envelope equal to 800 K. However, we do not adjust the 
inner radius since we do not need such rigorous adjustment for the purpose of galaxy modelling. 
As a result, our temperature is slightly higher (~ 830 K) due to the heating by the scattered 
light and dust emission. This mutual heating of dust is negligible when ri Mm < 100. On the 
other hand, in the case of p — 2, the effect of mutual heating is larger due to more pronounced 
concentration of the dust near the inner radius of the shell than the case of p = 0. In this 
case, the adjustment of the inner radius of dusty zone is indispensable to compare the results 
properly with the benchmark solutions. 

3. Sensitivity Analysis 

It is necessary to know how the resulting SEDs depend on the model parameters such 
as age t, optical depth ry, core radius ratio rj, size scaling factor R, and type of the extinction 
curve, in order to judge whether the results of SED fitting are unique or not. In this section, we 
describe how the radial temperature distributions and the SEDs depend on these parameters. 
The standard model is defined by a set of parameters t = 10 Myr, ry = 1.0, rj = 1000, R = 1.0, 
and the MW type extinction curve, as tabulated in Table 3. Unless otherwise stated, the 
standard model parameters are used throughout this section. 

3.1. The Effect of Age 

Figure 12a shows the radial temperature distribution T eq (r) of models for various ages 
of the system. Due to high concentration of stars, the energy density of the radiation field is 
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highest at the centre, and the dust temperature attains maximum there. As stellar populations 
aging, T eq (r) becomes systematically lower in a whole system. As will be mentioned below, this 
leads to a shift of dust emission toward longer wavelength. At the same time, the temperature 
gradient becomes flatter, since the effective wavelength of dust heating photons becomes longer. 
Such photons can travel easier to the outer regions and heat the dust there. However, this 
flattening of temperature distribution does not affect the SED, significantly. 

Figure 12b shows a behaviour of the SED as a function of age. The SED of dust-rich 
objects is composed of the stellar and dust emissions with the peaks in the UV-NIR and MIR- 
submm ranges, respectively. Both the stellar and dust emissions decrease as a system aging, 
because both the emitted UV photons and the absorbed ones by dust become less numerous. 
When a system is young, the dust emission dominates the SED, while it becomes less important 
for older ages. Simultaneously, the dust becomes cooler and the dust emission peak shifts from 
~40/im (for 10 Myr) to 100/im (for 10 Gyr). Note that neither the dust absorption nor the dust 
emission affects the NIR part of SED significantly, while the fraction of energy emitted in FIR 
decreases monotonically as a system aging. This indicates that the FIR luminosity normalized 
by the H-band luminosity, Lfir/Lh, is sensitive to age (see also Takagi et al. 1999), where 
Lpin is the FIR luminosity derived from 60 and 100 /im fluxes. 

3.2. The Effect of Optical Depth 

Figure 13a shows T eq {r) for different optical depths, Ty. T eq (0) is almost independent of 
the optical depth for ry < 40. When a mean free path of the UV photon is larger than the core 
radius of stellar distribution, almost all stars in the system can contribute to the heating of 
dust at the system centre. Note that for the homogeneous distribution of dust the optical depth 
within r C) s is 10 2 2 times smaller than that defined with r t in the adopted geometry. The SED 
of dust emission becomes broader with increasing Ty, because the temperature at the outer 
region drops sharply, while T eq (0) is kept constant. This is not true for the case of r\ < 10. In 
this case, T eq {r) is less sensitive to Ty, since strongly concentrated dust obscures the emission 
from the centre, and the energy density in the outer region is determined by stars in the outer 
region. 

The SEDs with various Ty are shown in Figure 13b. When ry increases, the stellar 
emission decreases due to the dust absorption, as a result of which the dust emission increases; 
thus shows a remarkable contrast to the age effect. The fraction of energy emitted in the FIR 
increases rapidly as a function of Ty when Ty <C 1, but this fraction is less sensitive for Ty > 1. 
Note that the peak wavelength of dust emission is less sensitive to Ty than to the age, since 
only the dust temperature in the outer region depends on Ty, while the age affects the dust 
temperature in a whole system. 

It is worth noting that the SED approaches the source function for large Ty (> 10). 
In particular, this becomes conspicuous in the UV region, since the extinction cross-section 
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increases with decreasing A. Thus, the SED in this region resembles the source function for 
large Ty. If stars distribute point-likely and are enshrouded by dust shell, the source function, 
\ in dusty medium is zero. Therefore, there should be a sharp cut-off in the SED, i.e. 
the emission in the UV-NIR spectral region of the population enshrouded by dust is virtually 
negligible (e.g. Figure 5a in Rowan- Robinson 1980). 

3.3. The Effect of Dust Geometry 

Figure 14a shows the effect of dust geometry r\ = r^ D /r cS on T eq {r). T eq (0) drops slightly 
in the case of rj = 1, since the mean free path of UV photons becomes comparable to r c s- On 
the other hand, T eq (r) rises as rj decreases in the outer region, since the optically thin region 
from which the UV photons can escape from a system enlarges. However, this effect is rather 
minor when compared with the effect of age and Ty. 

Figure 14b shows the effect of dust geometry on the SED. The peak of dust emission 
moves toward shorter wavelength for smaller rj, although T eq (r) is not very sensitive to rj. This 
is because the fraction of dust located near the centre where the dust is hottest increases with 
decreasing r]. Note that even if Ty is the same, the amount of dust changes when r) is different. 
When the dust mass with i] = 1000 is normalized to 1.0, the value of dust mass becomes 0.6, 
0.03 and 0.006 for rj = 100, 10 and 1, respectively. Since there is less amount of dust in a system 
with smaller rj, Lpm drops with decreasing rj. If the mass of dust is kept constant for different 
r/, then the L FIR is of the same order of magnitude. In the cases of r\ <C 100, the resulting 
SEDs are close to the intrinsic SED. The bump feature at 2175 A virtually disappears because 
the emission from the optically thin region dominates the resulting SED. Even if Ty ^> 1, the 
UV-NIR SED does not change significantly. Thus, the observed spread of starburst galaxies in 
two-colour diagrams can be explained by the stellar and dust distributions only with rj > 100 
(cf. Gordon et al. 1997). For r/ > 100, it is difficult to determine the value of r) precisely from 
the SED fitting alone, since the resulting SEDs are degenerate in geometry and optical depth. 
For instance, the SED of stellar emission with rj = 1000 and Ty = 1 can be reproduced by the 
model with rj = 100 and Ty slightly larger than 1. 

3.4- The Effect of System Size 

Figures 15a and 15b show T eq (r) and the SEDs for various scaling factors R, respectively. 
Note that the dust mass varies with R for constant Ty. The SED of absorbed stellar emission is 
independent of the physical size of a system, provided that Ty is kept constant. Since M* does 
not vary, T eq {r) rises as a whole with decreasing R due to the increase of energy density which 
leads to the shift of peak of dust emission toward shorter wavelength. The distribution of the 
source function remains the same, since both distributions of stars and dust depend on R in 
the same way. The scaling relation for absolute value of the source function is given by 
oce* /na e , where e* ocR~ 3 . Since Tocna e R, we get finally S^ 1 ' ocR~ 2 /t. Here A dependence 
is omitted for clarity. Since the resulting energy density of radiation field and have the 
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same dependence on R, the integral in the r.h.s. of the equation (5) is scaled by R~ 2 when the 
mutual heating of dust is negligible. That integral is proportional to T^ q +d , if oc \~ 6 . Thus, 
T eq (r) is roughly scaled by R~ 2 ^ 4+9 \ 

The effect of R on the peak wavelength of dust emission is comparable to that of the 
age. Once the age and ry are fixed, R could be determined by the peak wavelength of dust 
emission. This point is important since the results of SED fitting can be restricted by the size 
of the system derived from the observation. Thus, it can help to determine the dust geometry 
in the objects with Ty <C 1, like elliptical galaxies (see §4.3). 

3. 5. The Effect of Extinction Curve 

Figures 16a and 16b show T eq (r) and the SEDs for the MW, LMC and SMC extinction 
curve, respectively. In order to enlarge differences among the SEDs, Ty = 30 is adopted, instead 
of the standard value. The systematic difference in T eq {r) should be noted, due to the different 
strength of self-absorption indicated by the depth of the silicate absorption feature at 9.7 /im. 
The resulting shift of the peak of dust emission, however, is negligible in comparison with those 
due to age or R. Although the same Ty is adopted for the cases of all extinction curves, the 
resulting luminosity at V band is different because of the different amount of scattered light. 

3.6. The Fraction of Absorbed Energy 

We introduce the energy fraction absorbed by dust, AE/E. The resulting SEDs are 
integrated from 91.2 ran to 3 /im to calculate the amount of unabsorbed energy, E*, and then 
AE/E is calculated by using the total energy, E, and AE = E — E*. In Figure 17 and Table 
4, AE/E is presented for various ages, Ty, rj, and extinction curves. Note that AE/E is 
independent of R and is equal to the energy fraction emitted by dust. The dust distributed 
with larger rj absorbs the energy more efficiently. This is due to an increase of the effective 
amount of dust surrounding the concentrating stars. For the constant value of Ty and r], AE/E 
is always larger in younger system, since the absorbed UV photons are more abundant. The 
dependence on the extinction curve is negligible except for the cases of youngest age 10 Myr 
and Ty < 1. AE/E increases rapidly as a function of Ty when Ty < 1. This strong dependence 
becomes less pronounced for an old system. When Ty > 1, AE/E becomes more sensitive to 
the age rather than Ty, especially for young system. 

3. 7. Summary of Sensitivity Analysis 

We summarize the main results of our analysis: 1) the geometry of starburst galaxies 
should be described by 77 > 100, since the SED of absorbed stellar emission cannot be red 
enough to reproduce observed SEDs for 77 <C 100. 2) the UV-NIR SED is most sensitive to 
the stellar age and optical depth Ty when 77 > 100, 3) the FIR excess defined as Lfir/Lh is 
sensitive to the age but not to Ty when Ty > 1, which makes it useful for estimating the age of 
system independently of Ty, 4) the dust emission reaches maximum at the FIR-submm region, 
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and the peak wavelength of dust emission depends mainly on the age and R. From this study, 
we conclude that one can roughly estimate the stellar age from Lfir/Lh, tv from the shape 
of the UV-NIR SED, and R from the peak wavelength of dust emission The depth of silicate 
absorption feature at 9.7 /im is useful to constrain the type of extinction curve. By using the 
resulting R and the total stellar mass M* derived from the normalisation of SED, the size of 
stellar system is calculated, which can be compared with the observed size of the system to 
disentangle the degeneracy among R, ry and rj for the optically thin case (see §4.3). 

4. Comparison with Observations 

The best-fitting SED model to observed data has been sought by trial and error, using 
the results of sensitivity analysis given in the previous section. The values of x 2 are calculated 
for each fit. Although, in the SED fitting procedure, data at the MIR wavelength region are 
used to constrain 7y and to choose the plausible extinction curve, we exclude these data from 
the calculation of \ 2 ■ This is because the SED at this wavelength region is very sensitive to the 
adopted optical properties of silicate grains and PAHs, which are still a matter of debate, and 
therefore cannot be used to judge the goodness of the fit. For each galaxy, we show not only 
the best-fitting model, but also the near-fitting models, which are helpful to infer the accuracy 
of derived model parameters. 

As it is shown in §3.3, r] should be larger than 100 to reproduce the strongly attenuated 
SEDs of starburst galaxies, although it is difficult to constrain the precise value of r\ from 
the SED alone. Here, we adopt rj = 1000 for starburst galaxies, as a result of which the dust 
geometry is very close to the simple homogeneous distribution. The model parameters to be 
determined by observed SEDs are the age of stellar populations t, optical depth in ^-band, 7y, 
and the size scaling factor R. For an elliptical galaxy NGC2768, the SED from UV to NIR 
is reproduced by a galactic wind model, following Kodama & Arimoto (1997), and then the 
observed SED of dust emission is used to constrain ry, r\ and R. Thus, actually we fit SEDs 
of both starburst and elliptical galaxies by employing three free parameters. The resulting 
parameters for the best-fitting models and the derived quantities are tabulated in Tables 5 and 
6, respectively. 

4.1. Arp220 

Arp220 is the nearest ultraluminous infrared galaxy at a distance 77 Mpc away from 
the Milky Way (Soifer et al. 1987). Recent observations by ISO suggest that the source of 
dust heating is nuclear starburst rather than the AGN (Lutz et al. 1996; Lutz, Veilleux & 
Genzel 1999). Wynn- Williams & Becklin (1993) derived the "compactness" of MIR emission of 
Arp220. The compactness is defined by the ratio of fluxes observed in two different apertures, 
5". 7 and 100". They show that about 99% of MIR emission comes from the region within 5". 7 
aperture. Furthermore, Soifer et al. (1999) observed Arp220 with the Keck II telescope to find 
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that the 24.5 /im flux in the 4" diameter is consistent with the flux density at 25 /im obtained 
by IRAS in a rectangular beam, ~ 1' x 5'. Thus, MIR emission of Arp220 is concentrated in 
the central region. Here, we confront our models to the central area within 5" aperture, which 
corresponds to ~ 2.0 kpc. 

The photometric data are taken from Sanders et al. (1988; B, g, r, i), Carico et al. 
(1990; J, H, K, L), Smith, Aitken & Roche (1989; MIR), Klaas et al. (1997; MIR to FIR), 
NED (NASA/IPAC Extragalactic Database; IRAS bands) and Rigopoulou, Lowrence & Rowan- 
Robinson (1996; submm). Recently, Goldader et al. (2002) presented the far-UV (1457 A) and 
near-UV (2364 A) observations with HST/STIS. As noted by the authors, these photometric 
data are problematic because the large angular size of Arp220 is enough to fill the entire field 
of view of STIS, which results in the highly uncertain flux level of sky. Therefore, these data 
are not used for the SED fitting, although we show them in the related figures. 

In Figures 18a-d, we show 9 SED models in total, consisting of the best-fitting model 
and 8 near-fitting models. The parameters for these models are tabulated in Table 7 along with 
the value of reduced \ 2 ■ The best-fitting model is indicated by solid lines in all the four figures. 
It is found that the SMC extinction curve is required, in order to fit simultaneously the strong 
silicate absorption feature at 9.7 /im and the UV-NIR SED. Both the models with MW and 
LMC type extinction curve overestimate the MIR fluxes, when the UV-NIR SED is well fitted. 
Therefore, we use the SMC extinction curve for all the nine models. 

The best-fitting model is giving by t=30 Myr, ry = 40 and R = 0.5 with the SMC 
extinction curve, from which the bolometric luminosity L« = 1.4 10 12 L Q , the total stellar mass 
M* = 4.2 1O 1O M and the total dust mass Md = 4.8 10 8 M Q are derived. The resulting R = 0.5 
corresponds to r t = 3.2 kpc; thus, roughly consistent with the adopted aperture size (~ 2.0 
kpc). Note that, in order to compare the size of starburst region, we should be more precise for 
the definition of starburst region; for example, the effective radius observed at some wavelength 
should be compared with the effective radius of the model at the same wavelength. However, 
such a detailed comparison of the geometry is not our main purpose of this paper. 

The derived value of reduced \ 2 is relatively large even for the best-fitting model. In 
the case of Arp220 and also M82, a significant contribution to \ 2 comes from only few most 
deviating photometric data, e.g. at 3.7 /im for Arp220 and at 400 /im for M82. We would 
like to stress that the results of further fine-tuning of the SEDs parameters strongly depend 
on unknown systematic errors of observations performed with various telescopes, detectors 
and slightly different aperture sizes over a wide wavelength range. Indeed, genuine systematic 
uncertainties on data are expected. It is virtually impossible to obtain the exactly aperture- 
matched photometric data from the measurements performed with different aperture sizes and 
beam profiles for starburst regions whose surface brightness profile depends on the observed 
wavelengths (e.g. Scoville et al. 2000). Also, note that starburst regions are surrounded by 
underlying stellar populations with unclear boundaries. Furthermore, the absolute values of 
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X 2 depend on the accuracy of statistical errors which were estimated quite inhomogeneously. 
Therefore, we refer the value of \ 2 only to compare the goodness of fit among the reference 
models. 

The near-fitting models in Figure 18a have the same model parameters as the best-fitting 
model, but age. Similarly, the different ry and R are adopted for the near-fitting models in 
Figures 18b and 18c, respectively. The effects of age, ry and R are somewhat degenerate in 
this parameter range, especially due to large ry; i.e. 1) the variations of optical-NIR SED due 
to the age effects are very similar to those due to ry, 2) the average temperature of dust is 
changed by both the age and R effects. On the other hand, the UV luminosity is less sensitive 
in the range of age, t—15 - 60 Myr. This is because the bolometric luminosity is dominated by 
the UV luminosity, intrinsically. Note that, in this case, the bolometric luminosity is outlined 
by the FIR observations. 

In Figure 18d, we show the best fitting model together with the models of 10 Myr and 
300 Myr. Despite the large difference in age, the optical-NIR SEDs become similar by adjusting 
both ry and R. Therefore, the age of stellar populations is to be conservatively determined 
in the range of 10 - 300 Myr, when there are no constraints on ry and R from the other 
observations. 

Once the optical depth is determined independently of the SED, the age can be con- 
strained in the range of 15 - 60 Myr as shown in Figure 18a. In our model, ry is well restricted 
by the silicate absorption feature when ry ^ 10. For Arp220, the best-fitting is achieved with 
ry = 40 for t = 30 Myr. From the spectroscopic study with ISO, the extinction of ry ~ 45 
was derived for Arp220 (Sturm et al. 1996; Genzel et al. 1998). The observations of radio 
recombination lines result in the similar optical depth (Anantharamaiah et al. 2000). These 
observational results seem to be consistent with our result. However, it is misleading to compare 
these values directly, since they are derived by assuming the MW extinction curve. Therefore, 
we compare optical depth at the NIR region, where hydrogen recombination lines were observed 
by ISO to derive the extinction. In our model for the SMC type extinction curve, ry = 40 cor- 
responds to the optical depth at K band, tk = 2, while the observed value is tk = 4.5. This 
difference by a factor of 2 is not serious, since the observed optical depths are inferred from the 
emission lines of ionized gas, whose distribution could be different from that of dust. Calzetti, 
Kinney & Storchi-Bergmann (1994) suggest that the discrepancy by a factor of 2 is rather 
typical for starburst galaxies. Considering these results, we suggest that the constraint on age 
is rather tighter than the conservative case noted above. 

Among 9 models, only the oldest model (300 Myr) shown in Figure 18d is marginally 
consistent with far-UV data, but overestimates near-UV data. Although these data are less 
reliable as noted above, the observation by HST/STIS actually suggests very low fluxes at these 
wavelength, even possible uncertainty of the sky level is taken into account. Such a low flux at 
the far-UV may be due to either the selective extinction around massive stars or to unexpectedly 
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large anomaly of the extinction curve. When massive and hot stars tend to be embedded in the 
clouds whose optical depth is considerably larger than the average which we derive, the far-UV 
light dominated by these stars can easily be overestimated in our model. As for the extinction 
curve, it seems neither MW, LMC nor SMC extinction curves could be applicable to starburst 
regions, especially in the far-UV region in which very small grains play dominant role, since 
the highly energetic environments probably affect the size distribution of dust grains. In order 
to reconcile the best-fitting model with the observed far- and near-UV fluxes, the extinction 
should be extremely effective at these wavelengths by an order of magnitude comparing with 
the SMC extinction curve. These possibilities should be investigated carefully in the future 
with reliable observations. 

4.2. M82 

M82 is a nearby starburst galaxy at a distance of 3.3 Mpc. This galaxy shows remarkable 
activities including the disturbed dusty appearance and evidence for explosive star formation. 
This may be caused by an interaction with a neighbouring spiral galaxy M81 located at ~ 36 
kpc away in projection (Yun, Ho & Lo 1993; Ichikawa et al. 1995). The centre of this galaxy 
exhibits strong IR emission extending ~ 30" (Telesco & Harper 1980) which coincides with the 
effective radius of the bulge component (Ichikawa et al. 1995). Thus, we confront our models 
to the SED of central starburst region within ~ 30" (~ 0.5 kpc). 

The photometric data are taken from Johnson (1966; U, B, V, R, I), Ichikawa et al. 
(1995, J, H, K), Rice et al. (1988; IRAS bands), Telesco & Harper (1980; FIR), Jaffe, Becklin 
& Hildebrand (1984; submm), and Elias et al. (1978; submm). The spectrum from 2.4 jum to 
45 /im is obtained with the similar aperture size by the Short Wavelength Spectrometer (SWS) 
on board the ISO (Sturm et al. 2000). From the ISO-SWS data, we adopt the fluxes at 5, 16, 
20, 30, 40 fj,m for the calculation of x 2 ; at which the continuum emission dominates. 

In a similar manner to the case of Arp220, we show the best-fitting model and near- 
fitting models in Figure 19a-d. The parameters for these models are tabulated in Table 8, along 
with the resulting values of reduced \ 2 - The best fit is achieved with t=30 Myr, ry=8.5 and 
R=1.0. We choose the LMC extinction curve for this galaxy, judging from the depth of silicate 
absorption feature and the goodness of fit to the optical-NIR SED. Using this best-fitting 
model, we obtain L hol = 3.7 1O 1O L , M* = 1.1 1O 9 M and M D = 9.0 1O 6 M . The resulting R 
corresponds to r t =1.0 kpc, which is larger than the observed size of starburst region by a factor 
of ~4. 

The behaviours of SED when changing t, ry and R are almost the same as for the case 
of Arp220. The notable difference is found for the change of ry, due to the smaller value of ry 
for M82; i.e. the difference in ry affects only the wavelength region shorter than ~ 2fim. For a 
given age, the uncertainty of ry is at most Ary ~ 2 as shown in Figure 19b, while Ary ~ 10 
for Arp220. In the case of no restriction for both age and Ty, the conservative range of age is 
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10 - 100 Myr as shown in Figure 19d, which is tighter than in the case of Arp220 due to the 
smaller value of ry. 

The observed flux ratio of hydrogen recombination lines (Pa/3/Br7) indicates tk — 1-0 in 
the central region (Satyapal et al. 1995). We derive tk = 0.5 applying the LMC type extinction 
curve, and thus find again a factor 2 discrepancy between the optical depth derived from the 
stellar continuum and the line ratios. 

In our model, the optical depth at 9.7/im is 1.1 for Ty = 8.5 with the LMC extinction 
curve, and therefore the silicate absorption is important. On the contrary, Sturm et al. (2000) 
concluded that no evidence for the strong 9.7 /im absorption feature is found. They compared 
the spectrum of M82 with that of galactic reflection nebula NGC7023 which is expected to be 
free from the foreground extinction. There is no significant difference in the spectra, except 
for A > 12 /im. They succeeded to reproduce the spectrum of M82 up to A ~ 12 /im with the 
scaled spectrum of NGC7023 plus a power-law continuum which starts at 8.5 /xm. However, the 
solution obtained by their model may not be unique, since the effect of self-absorption could be 
compensated by the assumed flux of power-law continuum, which is a free parameter in their 
case, but constrained by the other regions of SED in our model. 

Thus, it is found that both the SEDs of M82 and Arp220 can be explained by the stellar 
populations of the same age. We note that, although the derived optical depth is different by 
a factor of ~5, Arp220 and M82 have the similar mass ratio of dust to stars, ~ 0.01. This also 
suggests that both Arp220 and M82 are at the similar evolutionary phase as a starburst galaxy. 
Therefore, the large difference between their SEDs is mainly caused by the difference in the 
geometry; i.e. the starburst region of Arp220 is much more compact than that of M82. This 
argument is consistent with the obtained values of R, which are restricted mainly by the peak 
wavelength of dust emission; i.e. the average temperature of dust, in other words, the energy 
density of radiation field and therefore the stellar density. 

4.3. HR10 

HR10 at z = 1.44 is one of famous EROs having R — K > 6 in the observed frame 
(Graham & Dey 1996; Cimatti et al. 1998). In general, extremely red colours are attributed 
either to very old stellar populations or to young ones heavily obscured by dust (Graham & 
Dey 1996). The latter case applies for HR10, because of conspicuous dust emission (Dey et 
al. 1999). The photometric data are all taken from Dey et al. (1999), except for ISO data 
(Elbaz et al. 2002). For this distant galaxy, a precise choice of aperture for starburst region is 
difficult, and therefore underlying stellar populations could contribute to the observed fluxes. 
Dey et al. (1999) note that optical SEDs (in rest frame) is similar to that of nearby ULIRGs 
rather than that of normal spiral galaxies which is expected for underlying stellar populations 
in gas-rich galaxies, like HR10. Therefore, we assume that observed fluxes are dominated by 
starburst stellar populations. 
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In Figure 20, we show the best-fitting model with the age of 300 Myr, and the near- 
fitting models with the age of 100 Myr and 1 Gyr. We choose the MW extinction curve which 
gives better fit to rest frame UV data. The nearest fits are sought by fitting data excluding the 
MIR data, which are also excluded in the calculation of reduced x 2 ( see Table 9), like in the 
cases of Arp220 and M82. 

For the best-fitting model with t = 300 Myr, ry = 7 and R—1.0, we derive M* = 
3.3 10 11 M o and L^oi = 1-4 1O 12 L . The absorbed energy corresponds to ~ 90% of the to- 
tal one. The dust mass Md = 1.2 1O 9 M is larger than that of the Arp220, despite a smaller 
Ty. The resulting R corresponds to r t =18.2 kpc, while the resolved size of this galaxy is ~ 0".9 
(Dey et al. 1999) corresponding to ~ kpc. 

For HR10, the acceptable range of age is systematically older than Arp220 and M82. 
The very red SED of HR10 at UV-optical cannot be reproduced by the model with t < 100 
Myr. The resulting mass ratio of dust to stars is smaller than that of Arp220 and M82 by a 
factor of 2. This is consistent with the old age of HR10 as a starburst; i.e. gas, the fuel of 
starburst, decreases with aging, while stellar mass increases as a result of star formation. 

The SED without MIR data can be reproduced with the wide range of starburst age 
(100 Myr - 1 Gyr) as shown in the figure. However, the oldest model is strongly rejected 
by MIR data. Note that, for most of starburst galaxies at high redshifts like those found by 
the submillimetre common-user bolometer array (SCUBA) survey, MIR data are not available, 
yet. Also, observations at rest frame NIR and/or at the emission peak of dust are extremely 
important for the tight constraint on the age of high-z starburst galaxies. 

44. NGC2768 

NGC2768, an elliptical galaxy (E6) at a distance of 21.5 Mpc, is one of a few ellipticals 
that were detected in submm wavelength (Wiklind & Henkel 1995). Photometric data are 
taken from Longo, Capaccioli & Ceriello (1991; far-UV), RC3 catalogue (de Vaucouleurs et 
al. 1991; UBV), Frogel et al. (1978; JHK), and Wiklind & Henkel (1995; submm). All data 
are reduced with the aperture (A) of total diameter (D) of the galaxy, i.e. A/D—l. Although 
Frogel et al. (1978) gave JHK data for A/D = 1, the aperture correction they adopted is not 
consistent with the optical photometry. Therefore, we have adjusted JHK data to the system 
of the optical data. In the far-UV bands, only the central part of galaxy was observed by the 
IUE satellite (Longo et al. 1991). Since the far-UV observations are performed in rectangular 
apertures 10" x 20", we convert them to the effective circular apertures. We have to note that 
such procedure is rather uncertain unless a true far-UV photometric profile of the galaxy is 
known. 

Since it is obvious that the stellar population in an elliptical galaxy is by no means the 
SSP, we calculate the unabsorbed SED of NGC2768 by using a galactic wind model that takes 
into account a realistic star formation history in a galaxy (Arimoto & Yoshii 1987; Kodama & 
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Arimoto 1997). Model parameters are the same as Kodama & Arimoto (1997) who explained 
the empirical colour- magnitude relation of ellipticals in Coma cluster as a sequence of mean 
stellar metallicity; i.e. t — 12 Gyr, tsF = 0.1 Gyr, ti n = 0.1 Gyr and x = 1.1, where t is a galactic 
age, tsF is a time scale of star formation, tj n is a time scale of gas infall, and x is a slope of the 
Salpeter-like IMF for the stellar mass interval 0.10 < m/M Q < 60. NGC2768 has M v = -21.8 
mag, from which we derive t gw = 0.24 Gyr for an epoch of galactic wind 1 . 

Figure 21 shows the best-fitting models for NGC2768. Best fits are achieved for 
(Ty,r),R) = (0.3,1,12) and (0.035,1000,0.7) as tabulated in Table 10 with the value of re- 
duced x 2 - Although we choose the MW extinction curve, the resulting SEDs are not sensitive 
to the extinction curve, since the amount of reddening required is very small. Obviously, the fit 
is not unique as the value of x 2 suggests, because the SEDs from optical to NIR are insensitive 
to dust in the case of Ty < 1. The degeneracy can be resolved if one takes into account a 
observed size of the galaxy. The effective radius of NGC2768 is r e = 76". 5 which corresponds 
to ~ 8 kpc (Peletier et al. 1990). Using the ratio r e /r c ~ 10 given by the adopted density 
distribution of stars, we obtain r e = 0.88 and 11.0 kpc for R = 1.0 and 12, respectively; each 
value of R is converted to the real size of galaxy by using M* = 2.2 1O 11 M for Lf, i = 9.0 1O 1O L . 
Thus, the model with R = 12 and rj = 1 fits better to the observation, giving = 2.9 1O 6 M . 

If all the dust in NGC2768 comes from stars via stellar wind, the amount of dust may be 
given as Mjj = 5ZM g , where 5 is dust-to- heavy element ratio, and M g is the mass of gas. The 
abundance of gas should be similar to the mean stellar metallicity of a galaxy. Our galactic 
wind model fitted for NGC2768 gives [Fe/H] = -0.09, or ~ O.8OZ ; thus Z = 0.016. After an 
epoch of galactic wind, the gas has been continuously supplied to the interstellar space from 
evolving stars, and the amount of gas should become as large as ~ 20% of total stellar mass 
(Arimoto 1989). This gives roughly M g ~ 0.2 x M* = 4.3 1O 1O M . Therefore, if the dust has 
always been kept in a galaxy, its amount should be as large as ~ 2.7 1O 8 M for S = 0.4. 
Then the average mass lose rate of dust during 12 Gyr is ~ 0.02 M yr _1 . The amount of dust 
we derive from the SED fitting is 2.9 1O 6 M , which is much smaller than the value expected 
and corresponds to the amount of dust accumulated during the last ~0.1 Gyr. Obviously, most 
of dust ever ejected from stars was lost to the extragalactic space and/or evaporated. 

5. Discussion 

5.1. Star formation history of starbursts 

In our model, the intrinsic stellar SEDs of the SSP are used. Since galaxies are composed 
of composite stellar populations, the elapsed time from the onset of starburst, or simply absolute 
age of starburst can be derived only if the star formation history is known. In the model of 
Efstathiou et al (2000), the star formation history was parameterized with a time scale for 

1 Kodama & Arimoto (1997) give different wind epoch, since they adopted slightly steeper IMF with x = 1.2 



23 



an exponential decay of star formation rate, which is estimated from the SED. For M82, they 
required two bursts of star formation with different time scales (2 Myr and 6 Myr). They 
derived the starburst ages of 26.5 and 16.5 Myr for each starburst event, respectively. On the 
other hand, in our model, the SED of M82 can be reproduced with a single starburst population; 
i.e. no need to assume multiple starburst events, although we derive the age very similar to 
that derived by Efstathiou et al. (2000). Note that Silva et al. (1998) assumed different star 
formation history of starbursts as well as the starburst age and applied to M82. This means 
that the derived star formation history depends on the model assumptions, strongly. Since, 
at least, these models can reproduce the observed SEDs of M82, we clearly need the other 
observational constraints, in order to estimate the star formation history of starbursts. The 
efficiency of star formation which is important clue to the star formation history can be inferred 
from the ratio of the indicator of star formation rate, such as the bolometric luminosity, to the 
total mass of gas in the starburst region (Sanders et al. 1988). We suggest that an evolutionary 
SED model in which the chemical evolution is taken into account, could be helpful to derive 
the absolute age of starburst, when the observed SED and gas mass are reproduced with the 
evolutionary SED model, simultaneously. 

5.2. Implications for high redshift starbursts 

Recently, various galaxy populations have been found as candidates of starburst galaxies 
at high redshifts. These populations can be classified by the wavelength ranges in which galaxy 
populations are defined; one class consists of galaxies whose characteristics of stellar emission 
are used as a criterion, such as Lyman-break galaxies and EROs, and for another class dust 
emission is used, like galaxies detected by ISO and/or SCUBA. It is important to investigate 
the relationship between these two classes of galaxy population, in order to reveal the nature 
of high-z starburst galaxies. Lyman-break galaxies are too faint at submm wavelengths to be 
detected with the currently available telescopes (e.g. Ouchi et al. 1999). EROs have been 
expected as optical counterparts for SCUBA galaxies, since their extreme red colours can be 
explained by the significant dust extinction. Although this prospect seems to be confirmed with 
the submm detection of a prototypical ERO HR10, Cimatti et al. (2002) however found that 
the star formation rates of star-forming EROs with R — K > 5 are an order of magnitude less 
than those of SCUBA galaxies. This result is consistent with the negative SCUBA detection 
for R-K EROs. Mohan et al. (2002) detected none of 27 EROs defined with R - K > 5 or 
I — K > 5, but one, in the submm wavelengths. Thus, it is not yet clear even whether there is 
a representative optical/NIR galaxy population for SCUBA galaxies. 

Great opportunities to study statistical properties of high-z starburst galaxies will be 
provided by the future infrared space missions, such as SIRTF-SWIRE (Franceschini et al. 
2002) and ASTRO-F (Nakagawa et al. 2001). In both missions, the wide range of wavelength 
from NIR to FIR will be covered, as a result of which both stellar and dust emission of high-z 
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starbursts can be observed. Therefore, these missions will be very important to investigate 
the relationship between the galaxy classes characterized by stellar and dust emissions. Our 
model is optimal for this analysis, since all regions of SED are physically interrelated, as it was 
demonstrated for the nearby starburst galaxies. Furthermore, starburst ages, stellar masses and 
optical depths could be derived for statistical sample of high- z starburst galaxies, as a result 
of model application to the data obtained with both missions. These results should provide 
important clues on the process of galaxy formation. 

6. Conclusions 

We have developed a code which solves the equation of radiative transfer by assuming 
spherical symmetry for arbitrary radial distributions of stars and dust. Isotropic multiple 
scattering is assumed and the self-absorption of re-emitted energy from dust is fully taken into 
account. The adopted dust model takes into account graphite and silicate grains, and PAHs. 
We successfully reproduce the MW, LMC and SMC type extinction curves as well as the spectra 
of cirrus emission in the MW. 

We adopt the geometry where stars and dust are distributed by the King model. The 
SEDs of starburst galaxies are reproduced well with the core radius ratio of dust to stars 
77 > 100; thus, dust distributes nearly uniformly and surrounds centrally concentrating stars. 
The UV-NIR SED is most sensitive to the age and optical depth ry. The FIR excess defined 
by Lfir/Lh depends mainly on age, and is rather insensitive to ry when ry > 1, which allows 
rough estimate of starburst age by using L FIR /L H . The peak wavelength of dust emission 
is in the FIR-submm region, and depends mainly on the age and size scaling factor R. We 
show that the extinction curve in starburst galaxies is close to that in the SMC (Arp220) or in 
the LMC (M82), judging from the MIR features when the other parts of SED are reproduced. 
Therefore, the non-MW type extinction curve should be considered, in order to model the SEDs 
of starburst galaxies. 

UV-submm SEDs of starburst galaxies, Arp220, M82 and HR10, can be reproduced 
consistently with the reddened starburst stellar populations. The derived age of stellar popu- 
lation in the starburst region of Arp220 is 30 Myr for an instantaneous starburst (SSP model). 
Although the same age is derived for M82, the obtained optical depth is significantly different 
(ry = 8.5 for M82, and 40 for Arp220). The restriction to the starburst age of M82 is tighter 
than that of Arp220, due to the smaller optical depth. For prototypical star-forming ERO 
HR10, we derive systematically older age (300 Myr) than Arp220 and M82, which is required 
to reproduce the extremely red colour at optical-NIR region. This difference in evolutionary 
phase starburst 

is also suggested by the smaller mass ratio of dust to stars in HR10 than that of Arp220 
and M82. 

The optically thin limit is studied in case of the elliptical galaxy NGC2768, where ry, 
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galaxy size, and r) are degenerate. Determination of the galactic size from observation is required 
to disentangle this degeneracy. From the observed size and SED of this galaxy, the distribution 
of dust is found to be similar to that of stars. 
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Tabic 1 . Parameters of PAH features 



\\pm] i>o[cm 1 ] T[cm x ] A[\Q 25 cm 3 ) 

3.3 3030 23J5 1A N H 

6.2 1610 36.5 0.7 N c 

7.7 1300 59.0 2.0 N c 

8.6 1160 34.5 1.2 Nh 

11.3.... 885 17.0 14.1 N H 

12.7.... 787 32.5 11.1 N H 



Note. Nc and Nh are the number of carbon and hydrogen atoms in a PAH molecule, respectively. 



Table 2. The relative fractions of dust particles 





M D /M H 


Graphite 


Silicate 


PAH 






BG 


VSG 


BG 


VSG 




MW 


8.0 10" 3 


0.37 


0.08 


0.43 


0.07 


0.05 


LMC 


3.6 10™ 3 


0.10 


0.02 


0.74 


0.13 


0.01 


SMC 


1.9 10~ 3 


0.04 


0.009 


0.808 


0.138 


0.005 



Note. Mo/Mh is the mass ratio of dust to hydrogen. 



Table 3. Parameters of the standard model 



Age 


T V 


R 


■n 


Extinction Curve 


10 Myr 


1.0 


1.0 


1000 


MW 
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Table 4. Fraction (%) of the absorbed energy 



MW Extinction Curve 









10 Myr 








1 Gyr 








10 Gyr 




TV 


It 


lot 


loot 


iooot 


1 


10 


100 


1000 


1 


10 


100 


1000 


0.5 .... 


6 


26 


45 


49 


4 


11 


18 


20 


4 


9 


14 


15 


1.0 .... 


10 


39 


GG 


71 


5 


18 


31 


34 


5 


14 


24 


2G 


2.0 .... 


15 


52 


83 


87 


8 


29 


49 


53 


7 


24 


40 


44 


5.0 .... 


23 


GG 


93 


9G 


14 


4G 


74 


78 


12 


40 


GG 


71 


10.0.... 


30 


74 


97 


98 


19 


58 


87 


90 


17 


53 


82 


8G 



SMC Extinction Curve 









10 Myr 








1 Gyr 








10 Gyr 




TV 


1 


10 


100 


1000 


1 


10 


100 


1000 


1 


10 


100 


1000 


0.5 .... 


8 


32 


55 


59 


4 


11 


19 


20 


4 


9 


14 


15 


1.0 .... 


13 


45 


72 


77 


G 


19 


32 


34 


5 


14 


24 


2G 


2.0 .... 


18 


5G 


84 


88 


8 


29 


49 


53 


7 


23 


40 


43 


5.0 .... 


2G 


G9 


94 


9G 


14 


4G 


73 


77 


12 


40 


G5 


70 


10.0.... 


33 


76 


97 


98 


19 


57 


86 


89 


17 


52 


81 


85 



f Core radius ratio, rj 

Table 5. Parameters of the best-fitting models 



Name 


t [Gyr] 


TV 


R 


V 


Extinction Curve 


Reduced % 2 


Arp220 


0.03 


40 


0.5 


1000 


SMC 


1.69 


M82 


0.03 


8.5 


1.0 


1000 


LMC 


1.85 


HR10 


0.3 


7.0 


1.0 


1000 


MW 


1.64 


NGC2768 


12° 


0.3 


12. 


1 


MW 


1.67 



a assumed parameters 



Table 6. Derived galaxy luminosities and masses 



Name 


M* 


Lboi 


M D 




[lO n M ] 


[10 n £ Q ] 


[1O 7 M ] 


Arp220 


0.42 


14. 


48. 


M82 


0.011 


0.37 


0.90 


HR10 


3.3 


14. 


120. 


NGC2768 


2.2 


0.9 


0.29 
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Table 7. Model parameters for Arp220 



Name 


t [Myr] 


TV 


R 


V 


Extinction curve 


Reduced 


model 1 


30 


40 


0.5 


1000 


SMC 


1.69 


model 2 


15 


40 


0.5 


1000 


SMC 


2.94 


model 3 


60 


40 


0.5 


1000 


SMC 


5.37 


model 4 


30 


60 


0.5 


1000 


SMC 


2.34 


model 5 


30 


30 


0.5 


1000 


SMC 


4.63 


model 6 


30 


40 


0.3 


1000 


SMC 


3.62 


model 7 


30 


40 


1.0 


1000 


SMC 


2.74 


model 8 


300 


70 


0.1 


1000 


SMC 


3.06 


model 9 


10 


30 


1.0 


1000 


SMC 


3.98 



Table 8. Model parameters for M82 



Name t [Myr] ry R i] Extinction curve Reduced x 2 



model 1 


30 


8.5 


1.0 


1000 


LMC 


1.85 


model 2 


15 


8.5 


1.0 


1000 


LMC 


3.35 


model 3 


60 


8.5 


1.0 


1000 


LMC 


8.83 


model 4 


30 


10. 


1.0 


1000 


LMC 


2.51 


model 5 


30 


7.0 


1.0 


1000 


LMC 


2.64 


model 6 


30 


8.5 


1.5 


1000 


LMC 


6.95 


model 7 


30 


8.5 


0.5 


1000 


LMC 


3.75 


model 8 


10 


8.0 


0.1 


1000 


LMC 


2.72 


model 9 


100 


10. 


1.0 


1000 


LMC 


3.72 



Table 9. Model parameters for HR10 



Name 


t [Myr] 


TV 


R 


V 


Extinction curve 


Reduced x 2 


model 1 


300 


7.0 


1.0 


1000 


MW 


1.67 


model 2 


100 


5.0 


3.0 


1000 


MW 


3.59 


model 3 


1000 


20. 


0.1 


1000 


MW 


2.08 
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Table 10. Model parameters for NGC2768 



Name 


t [Gyr] 


TV 


R 


f] 


Extinction curve 


Reduced \ 2 


model 1 


12° 


0.3 


12 


1 


MW" 


1.64 


model 2 


12 a 


0.035 


0.7 


1000 


MW a 


1.69 



a assumed parameters 




0.1 1.0 10.0 

X [/-tm] 

Fig. 1. The spectral energy distributions (SEDs) of the simple stellar populations (SSP) (Kodama & 
Arimoto 1997) are given for ages: 10 Myr (top), 100 Myr, 1 Gyr, and 10 Gyr (bottom). The luminosity 
is given for the SSP with the initial mass (i.e. mass at zero age) of 1 M . 
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Fig. 2. The MW, LMC and SMC type extinction curves per hydrogen atom are shown in panel a), b) 
and c), respectively. The dot-dashed curve is for graphite, dotted for silicate, dashed for PAH, and solid 
for total extinction. The average extinction curves of each type are taken from Pei (1992) and Rieke & 
Lebofsky (1985), and are depicted by open circles. 
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Fig. 3. The spectra of cirrus emission in the MW. The line symbols are the same as in Figure 2. COBE 
observations are represented by open circles (Dwek et al. 1997). The other data are taken from Desert et 
al. (1990) and references therein. 




Fig. 4. The normalized dust density profiles given for various rj values. The horizontal axis is normalized 

by n (=n,s)- 
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Fig. 5. Scheme of dusty galaxy. The intensity is calculated and stored at the positions of filled circles. 
bj indicates the impact parameter for the j-th ray. The angular integration at ry, where j' = j + 1, 



performed by using the rays depicted by dashed lines and that with j 
defined by v and a radius vector. 



1. The cosine grid fj,j at Tj 



is 



IS 
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2 4 6 8 10 12 



Fig. 8. The fraction of light which escapes from a dust-rich galaxy. Stars and dust are assumed to 
distribute homogeneously in this case. We adopt u> = 0.6. The thick lines indicate the results of present 
model, while the thin lines are taken from Witt et al. (1992). Solid, dot-dashed, and dotted lines show 
the total light, the direct stellar light, and the scattered light, respectively. The analytical solution for the 
direct stellar light is shown by open circles. 




2 4 6 8 10 12 14 



Fig. 9. The ratio of scattered light to the direct stellar light as a function of ry for different number of 
the scattering terms n = 1 (bottom) to n = 10 (top) with An = 1. Note that lines for n > 7 are nearly 
identical. The results of Witt et al. (1992) are plotted as open circles. 
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Fig. 10. The effects of self-absorption. (a) The temperature distributions of graphite BG with 
a = 0.15 fxm. The solid and dashed lines indicate T eq (r) with and without the self-absorption, respec- 
tively, (b) The SEDs. The solid and dashed lines indicate \L\ with and without the self-absorption, 
respectively. The dotted line gives the intrinsic SED. The adopted parameters are rj = 1000, ry = 30.0, 
t = 10 Myr, and the SMC extinction curve. 
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3.0 




1 10 100 1000 

r/r, 




J. [urn] 

Fig. 11. A comparison with output of DUSTY code by Ivezic et al. (1997). Solid and dotted lines indicate 
the results of present model and those of Ivezic et al., respectively, (a) The temperature distribution of 
dust. The radius is given in units of a radius of central dust free cavity, r\. (b) The SEDs normalized by 
an integrated flux, F. 

1000F j p I 




r At A [fj.m] 

Fig. 12. The effects of age. Parameters, except for the age, are the same as for the standard model (see 
text), (a) The temperature distribution of graphite BG with a = 0.15 fim for 10 Myr, 100 Myr, 1 Gyr, 
and 10 Gyr old. (b) The same as (a), but for the SEDs. Dotted lines give the intrinsic SEDs. 
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Fig. 13. The effects of optical depth. Parameters, except for Ty , are the same as for the standard model 
(see text), (a) The temperature distribution of graphite BG with a = 0.15 /im for ry = 0.1 (top), 1.0, 3.0, 
10 and 40 (bottom), (b) The same as (a), but for the SEDs. The dotted line shows the intrinsic SED. 




Fig. 14. The effects of dust geometry. Parameters, except for rj, are the same as for the standard model 
(see text), (a) The temperature distribution of graphite BG with a = 0.15 /im for rj = 1 (dot-dashed line), 
10 (dashed line), 100 (dotted line), and 1000 (solid line), (b) The same as (a), but for the SEDs. The 
dotted line shows the intrinsic SED. Figures indicate r\ values which increase from top to bottom for the 
stellar emission, but from bottom to top in FIR-submm range. 
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Fig. 15. The effects of system size. Parameters, except for R, are the same as for the standard model 
(see text), (a) The temperature distribution of graphite BG with a = 0.15 /im for R = 0.2 (top), 1.0, 5.0 
and 10.0 (bottom), (b) The same as (a), but for the SEDs. The dotted line indicates the intrinsic SED. 




Fig. 16. The effects of extinction curve. Parameters, except for ry = 30 and the extinction curve, are 
the same as for the standard model (see text), (a) The temperature distribution of graphite BG with 
a = 0.15 /xm for the MW (solid line), the LMC (dashed line), and the SMC (dot-dashed line) extinction 
curves, (b) The same as (a), but for the SEDs. The dotted line indicates the intrinsic SED. 
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2 4 6 8 1 

Fig. 17. The absorbed energy fraction as a function of optical depth. Solid and dotted lines correspond 
to 10 Myr and 10 Gyr, respectively. The four lines for each age correspond to geometry with rj = 1 
(bottom), 10, 100, and 1000 (top). 
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Fig. 18. The results of SED fitting for Arp220. The SMC type extinction curve and 77 = 1000 arc adopted 
for all models. The other parameters for each model are denoted in the figure. We use solid circles for 
data which are used to calculate the value of \ 2 , otherwise open circles are adopted. The upper solid line 
at UV-NIR region is the intrinsic stellar SED with the age of 30 Myr. In all panels, the best-fitting model 
is indicated by solid lines. Horizontal axis is the rest frame wavelength, (a) The SED models with the 
same parameters, but age, as that of best-fitting model indicated by solid line, (b) The same as (a), but 
ry is varied, instead of age. (c) The same as (a), but R is varied, instead of age. (d) The fitting model 
with the fixed age (10 and 300 Myr), shown with the best-fitting model. 
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up n nr 




0.1 1.0 10.0 100.0 1000.0 0.1 1.0 10.0 100.0 1000.0 

\ [fj.m] \ [M m ] 



up n up 




X [/lim] \ [M m ] 

Fig. 19. The same as Figure 19, but for M82. LMC type extinction curve and r\ = 1000 arc adopted for 
all models. Thick dashed line is for the ISO-SWS spectrum. The upper solid line at UV-NIR region is the 
intrinsic stellar SED with the age of 30 Myr. 




Fig. 20. The same as Figure 18d, but for the high redshift galaxy, HR10. ISO data are indicated by 
open circles. The MW extinction curve and r\ = 1000 are adopted for all models. The upper solid line in 
UV-NIR is the intrinsic stellar SED with t = 300 Myr. 
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Fig. 21. The results of SED fitting for the elliptical galaxy, NGC2768. The MW type extinction curve 
and t = 12 Gyr are adopted for all models. The stellar emission (12 Gyr) in FIR is indicated by lower 
solid line. 
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